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Abstract: 

A  new  computer  program  for  accurate  calculation  of  acoustic  ray  paths  through  a  range-varying 
ocean  sound  channel  has  been  written.  It  is  based  on  creating  a  model  of  the  speed  of  sound  in  the 
ocean,  consistent  with  input  data,  that  produces  the  smoothest  possible  wavefronts.  This  scheme 
eliminates  “false  caustics”  from  the  wavefront.  It  may  be  useful  in  calculating  an  approximate 
solution  to  the  full  wave  equation  at  megameter  ranges. 
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1.  Introduction 


The  Ray  program  is  a  part  of  an  ongoing  effort  by  John  Spies berger’s  research  group  to  create 
a  fully  automated  system  for  performing  basin  scale  ocean  acoustic  tomography. 

One  key  part  of  such  a  system  is  forward  modeling  of  multipath,  which,  given  a  model  ocean 
sound-speed  field  and  bathymetry,  predicts  when  sound  emitted  from  a  source  will  arrive  at  a 
receiver  and  where  the  sound  will  travel.  There  are  currently  several  ways  to  solve  the  forward 
problem  including  normal  modes,  the  parabolic  equation  and  geometric  ray  tracing.  Some  theoret¬ 
ical  investigations  into  forward  modeling  have  led  us  to  the  concept  ofeigentubes  which  are  the  full 
wave  generalization  of  eigenrays  [Bowlin,  1991).  These  considerations  lead  to  the  conclusion  that 
the  wavefront  generated  by  a  geometrical  ray  trace  is  useful  in  solving  the  forward  problem  at  fre¬ 
quencies  much  lower  than  would  be  expected  by  a  simple  analysis  of  the  high  frequency  assumptions 
implicit  in  the  geometrical  acoustics  approximation. 

The  idea  we  use  is  that  for  a  given  physical  situation  (environment  and  source  placement)  the 
geometrical  wavefront  will  be  independent  of  the  frequency  and  bandwidth  of  the  source.  Consider 
a  collection  of  geometric  rays  all  starting  from  a  single  source  location  wit  h  different  starting  angles 
and  all  ending  at  the  receiver  range.  The  depths,  arrival  times,  and  angles  of  these  rays  (  at  the 
receiver  range)  as  a  function  of  launch  angle  are  what  we  call  a  wavefront.  The  actual  acoustic  field 
for  a  given  source  frequency  and  bandwidth  can  then  be  calculated  from  the  geometrical  wavefront 
[Buchal  and  Keller,  I960).  We  envision  the  geometrical  rays  as  a  (frequency  independent)  skeleton 
to  which  a  flesh  is  added  with  (frequency  dependent)  diffraction  effects. 

In  addition  to  offering  a  framework  for  calculating  the  complete  acoustic  field,  the  geometrical 
optics  or  ray  model  offers  accurate  and  numerically  efficient  determination  of  distinctly  separated 
multipath  signals.  The  multipath  calculation  is  equivalent  to  finding  the  spatial  structure  of  a 
wavefront  formed  from  an  impulse  source,  a  structure  considered  since  the  early  stages  of  acoustic 
tomography  [Brown  et  al.,  1980).  This  structure  was  recently  analyzed  with  an  experiment  [ Duda 
et  al.,  1992). 

None  of  the  existing  ray  trace  codes  of  which  we  are  aware  are  suitable  for  the  full- wave  numer¬ 
ical  extension  at  the  megameter  ranges  of  basin  scale  tomography.  Our  primary  considerations  for 
evaluating  ray  trace  codes  have  been  accuracy  and  speed.  To  add  a  diffracted  flesh  to  the  geometric 
bones,  we  must  also  include  “smoothness  of  the  wavefront”  as  a  primary  consideration. 

Most  fast  ray  tracing  codes  for  ocean  acoustics  approximate  the  sound  speed  of  the  ocean  as 
piecewise  linear.  The  discontinuities  in  the  first  derivative  of  these  piecewise  linear  approximations 
cause  “fjse  caustics”  or  more  precisely  “non-turning  point  caustics”  (NTP  caustics).  One  way 
to  see  this  is  to  look  at  the  group  velocity  as  a  function  of  launch  angle  in  a  range  independent 
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environment.  Denote 


where  i>  is  the  launch  angle  and  va  is  the  “one  loop  group  velocity”  defined  as  the  range  of  a  single 
loop  of  a  ray  divided  by  the  time  it  takes  the  ray  to  traverse  the  loop.  A  little  bit  of  calculus  shows 
that  5  depends  upon  integrals  along  the  ray  of  the  form 


dz 
sin  # 


where  a,  and  (3 ,  are  bounded,  smoothly  varying  functions  and  d  is  the  angle  of  the  ray  with  respect 
to  the  horizontal.  For  a  piecewise  linear  sound  speed,  ^/(c7)2  becomes  a  delta  function  while 
1/sintf  has  an  integrable  infinity  (with  respect  to  dz)  at  the  turning  points  of  the  rays.  When  the 
turning  point  of  a  ray  approaches  a  knot  in  a  piecewise  linear  ocean  then  5  can  become  arbitrarily 
large.  This  causes  discontinuous  jumps  in  the  wavefront  which  would  not  occur  if  the  second 
derivative  of  the  sound  speed  were  bounded.  Caustics  appear  in  a  wavefront  at  extrema  of  the 
function  of  depth  at  the  receiver  range  vs.  launch  angle,  where 


This  condition  is  satisfied  whenever  the  wavefront  contains  a  ray  at  a  turning  point.  These  we  call 
turning  point  caustics.  Wherever  there  is  a  jump  in  the  wavefront  due  to  a  very  large  value  of  S 
then  a  caustic  can  also  arise.  These  caustics  are  not  restricted  to  lie  on  a  turning  point  of  a  ray  so 
we  call  them  NTP  caustics. 

Measurements  of  acoustic  pulses  at  a  few  hundred  Hz  and  at  distances  of  one  to  three  thousand 
kilometers  show  that  the  wavefront  is  simpler  than  predicted  by  ray  traces  which  model  the  sound 
speed  as  piecewise  linear  [Duda  et  al.,  1992;  Sparrock,  1990;  Spiesberger  and  Metzger ,  1991].  A 
reasonable  requirement  of  a  ray  calculation  (or  ray-tracing)  code  is  the  reproduction  of  a  stable 
averaged  wavefront  structure  when  a  smooth  averaged  ocean  is  considered.  The  program  Ray  is 
designed  to  produce  continuous  wavefronts  at  megameter  ranges.  The  overall  philosophy  has  been 
to  find  a  simple  model  of  the  environment,  consistent  with  the  input  environmental  data  (sound 
speeds  and  depths),  that  produces  the  smoothest  wavefronts.  This  strategy  eliminates  many  but 
not  all  of  the  NTP  caustics.  The  NTP  caustics  that  remain  are  due  the  the  sound-speed  structure 
of  the  environment.  These  irreducible  NTP  caustics  extend  smoothly  over  a  finite  fraction  of  the 
wavefront.  This  means  that  the  acoustic  field  in  the  shadow  zones  associated  with  these  caustics 
can  be  found  analytically  from  the  geometrical  wavefront. 

Although  the  idea  of  smoothing  sound-speed  profiles  in  order  to  produce  smooth  wavefronts 
has  been  around  for  a  long  time  ( Pederson ,  1961],  we  have  implemented  an  automated  version  of 
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the  smoothing  that  is  integrated  with  efficient  numerical  techniques  that  enable  us  to  trace  rays 
quickly  through  a  realistic  model  of  the  ocean. 

Mathematical  algorithms  of  Ray  are  described  in  section  2.  Sections  3,  4  and  5,  respectively, 
describe  the  input,  the  operation,  and  the  output  of  the  program.  Section  6  briefly  describes  the 
performance  of  Ray.  Section  7  is  a  summary. 
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2.  Computational  Algorithms 


Equations  of  Motion  and  the  Spherical  Earth  Correction 


The  equations  of  motion  for  a  ray  travelling  through  the  ocean  cam  be  cast  in  Cartesian  coor¬ 
dinates  as  follows, 

dr  c  c 

dz  . 

—  =  tan  9 
dr 

df  _  sec0 
dr  c 

where  9  is  the  angle  of  the  ray  with  respect  to  the  horizontal  r  axis,  and  z  is  the  vertical  coordinate. 
These  equations  are  derived  from  Fermat’s  principle  of  least  time  in  appendix  A. 

For  long  range  ocean  acoustics,  the  curvature  of  the  Earth’s  surface  makes  non-Cartesian 
coordinates  more  suitable  for  ray  tracing.  Let  new  i  axes  lie  along  radii  passing  through  the  center 
of  the  Earth  with  z  —  0  at  sea  level  and  z  =  Rc  at  the  earth’s  center,  where  Rt  is  the  radius  of 
the  earth,  and  let  the  new  f  be  the  range  measured  along  a  circular  arc  at  sea  level.  Three  things 
happen  when  the  equations  of  motion  are  translated  into  this  new  coordinate  system.  The  first  is 
a  trivial  change  in  the  sign  of  z  and  9.  The  second  effect  is  from  the  conversion  from  dr  to  dr. 
Imagine  a  fish  that  stays  at  a  depth  i,  directly  below  a  moving  boat.  The  fish  will  travel  a  shorter 
distance  than  the  boat  by  a  factor  of  ft  =  dr/df  =  ( Re  -  i)/Re.  The  third  effect  is  a  rotation  of 
the  coordinate  system  by  df/Re  radians  as  we  take  a  step  of  size  dr. 

The  new  equations  of  motion  which  include  geometrical  effects  due  to  a  spherical  earth  are 


d6_  die 
df~  e  c 

=  fe  tan  9 
dr  Je 

dt  _  fe  sec  6 
dr  c 


These  are  the  equations  that  Ray  integrates  subject  to  the  modifications  and  approximations 
discussed  below. 


Smoothing  the  sound  speed 

The  interpolation  of  sound  speed  as  a  function  of  depth  is  the  most  important  algorithm  in 
Ray.  For  accurate  calculation  of  ray  paths  and  ray  travel-times,  sound  speed  is  required  at  arbitrary 
locations.  A  procedure  is  required  which  takes  as  input  a  set  of  sound  speeds  defined  on  a  discrete 
grid  of  depths,  the  most  prevalent  form  of  input,  and  produces  a  function  describing  the  sound 
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speed  at  any  depth.  The  sound-speed  field  should  be  continuous  with  a  bounded  second  derivative, 
and  should  he  a  realistic  approximation  of  the  ocean. 

Begin  with  a  range  independent  sound-speed  profile.  Let  the  input  depth  grid  and  sound 
speeds  be  labeled  z,  and  c,  for  1  <  i  <  N  such  that 


c(Zi)  =  c,. 

Each  depth  z,  is  called  a  knot  for  reasons  that  become  apparent  below.  The  continuous  piecewise 
linear  approximation  to  these  points  is  denoted  c^l(z).  It  will  be 

c(1)(*)  =  c.  +  0i{z  ~  Zi),  Zi<Z<  zl+1 


with  s  (c,+i  -  e,)/(zl+i  -  z,).  This  is  a  continuous  function.  The  first  derivative  is  piecewise 
constant  of  the  form 


dc(*) 

dz 


=  A. 


Z,  <  Z  <  Z,+1. 


The  second  derivative  consists  of  a  sum  of  delta  functions, 


dV1* 

dz2 


N-X 

=  a'6(z  ~ 
i=2 


where  a,  s  /?,  -  To  smooth  out  the  the  sharp  corners  at  the  knots,  and  to  reduce  the 

magnitude  of  the  second  derivative,  one  can  convolve  this  continuous  piecewise  linear  model  with 
a  normalized  symmetric  tophat  function  defined  by 


9(w;z) 


M  <  v, 
\z\  >  w 


The  result  of  this  operation, 

c(3)(u>;z)s  j  c^\z  -  z')g{w\z')dz' 


is  a  continuous  piecewise  parabolic  function  with  a  continuous  piecewise  first  derivative  and  a  finite 
piecewise  constant  (but  discontinuous)  second  derivative.  If  the  width,  tv,  o*  g( )  is  less  than  the 
spacing  between  knots,  then  will  be  parabolic  for  z  within  w  of  a  knot  and  it  will  be  linear  and 
equal  to  for  all  z  that  are  not  within  w  of  a  knot.  In  the  parabolic  regions  the  second  derivative 
of  c<J>  is  constant,  and  in  the  linear  regions  it  is  zero. 

There  is  a  tradeoff  in  determining  w,  the  width  of  g(w;  z).  We  are  forced  to  use  an  interpolation 
scheme  because  of  our  ignorance  of  the  actual  c(z)  away  from  the  points  z,  where  c  is  measured. 
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A  w  is  sought  which  gives  the  smoothest  wavefront  possible,  and  still  remains  consistent  with  the 
original  set  of  measured  values,  {c,}.  The  wider  the  u>,  then  the  smaller  the  second  derivative  of 
and  the  smoother  the  wavefront.  For  equally  spaced  z,  one  choice  is  to  set  vi  =  (2,+i  -  z,)/ 2, 
which  will  eliminate  all  of  the  linear  sections  of  c*2*  where  the  second  derivative  is  zero  and  tend 
to  minimize  the  second  derivative  in  the  parabolic  sections. 

A  complete  set  of  data  that  is  available  for  determining  global  ocean  sound  sound  is  the  Levitus 
data  base  of  temperature  and  salinity  [Levitus,  1982],  which  does  not  have  equally  spaced  z,.  The 
above  method  of  smoothing  can  be  generalized  to  let  the  parameter  w  vary  with  depth.  Define  a 
set  of  widths  u>,,  such  that  the  width  of  g(z)  is  u>,  when  z  =  2,.  If  the  constraints 

2,  -  tu,  >  2,_1 
2.  +  «>,  <  2I  +  1 

are  imposed  on  the  set  of  widths,  then  it  is  always  possible  to  use  these  widths  to  construct  a 
continuous  piecewise  parabolic  function  (with  continuous  first  derivative)  analogous  to  c*2),  which 
will  be  equal  to  c(2)(u>,;2)  for  2  such  that  max(2,_i  -f  tc,_i, r,  —  w,)  <  2  <  min(2l+i  -te,+i,Zi  +  u>,). 
This  is  not  proven  here,  instead  expressions  arc  given  for  the  resulting  smoothed  function,  denoted 
by  c<3). 

The  functional  form  of  c<3)  will  depend  on  the  values  2,,  and  w,.  There  are  three  cases.  Case 
1  reproduces  the  linear  segments  of  c*1)  far  away  from  knots.  It  occurs  for  all  z  such  that 

2,-1  +  ui,_j  <  2  <  z,  -  w, 

then 

C<3>(2)  =  C,  -1-  (3,(z-  2,). 

Case  2  reproduces  the  parabolic  sections  of  c<2l(ut,;  2)  near  a  knot.  It  occurs  for  all  2  such  that 


max(z,_i  +  u/,^,2,  -  Wi)  <  .  <  min(2l+i  -  u\+i,2,  +  w,) 


then 


C<3>(2)  =  C,  + 


+  yh±L±Jhl(z  r.)  +  _2i_r,  _  ,.12 

2  (Z  Z'>+4w,{Z  Z'}- 


Case  3  is  the  smooth  interpolation  of  two  “overlapping”  sections  of  case  2.  It  occurs  for  all  2  such 
that 


2,+i  -  u»i+i  <  2  <  2,  +  u;,. 

Note  that  every  section  of  case  3  is  sandwiched  between  two  case  2  sections.  Let  c^3l( z )  defined  on 
each  of  these  section  be  0/(2)  and  cr(z).  Then 

c(3)(z)  =  c,(z,)  +  7,(2  -  2,)  +  — ~~(z  -  2, )2 

2(2r  -  2,) 
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where 


z!  =  ■*{+!  -  Wi+1 


=  2i  -  U-. 


7/  = 


7r  - 


dc; 


An  algorithm  to  automatically  generate  the  optimized  widths  for  any  depth  grid  has  not  been 
developed.  Ray’s  default  set  of  {w;}  has  been  designed  to  work  with  sound-speed  profiles  on  the 
depth  grid  of  the  Levitus  data  base. 


Removing  the  bios  created  by  smoothing 


There  is  a  serious  side  effect  of  the  smoothing  process  described  abo.e.  The  resulting  smooth 
function  does  not  pass  through  the  original  data  set, 


ew(*<)  =  c,  +  -~ 


If  the  widths  of  the  smoothing  regions  w,  were  much  smaller  than  the  spacing  of  the  z,,  then  this 
would  be  a  minor  problem  since  it  would  be  a  small  correction  over  a  small  region  of  z.  We  want 
to  make  the  widths  as  large  as  possible,  so  it  is  necessary  to  fix  the  discrepancy  between  the  input 
data  set  and  the  smoothed  profile  in  order  to  minimize  the  bias  in  the  travel  times  of  calculated 
rays. 

It  is  always  possible  to  find  a  new  set  of  sound  speeds  {c, }  which,  when  smoothed,  provide  a 
function  that  goes  through  the  original  data  points.  The  new  set  of  sound  speeds  can  be  determined 
by  solving  the  following  set  of  linear  equations  for  {<?,} 

„  2  i  W<  /^l  +  l  ~  ~  ‘V  i 

C|  -  C,f  —  !  - - - - - 

4  \  A,  A,_i 


with  A,  =  z,+i  -  2| .  Ray  does  not  solve  this  set  of  equations  directly.  Instead,  it  provides  a  method 
for  obtaining  an  iterated  solution  of  the  form 


4A,A,_ic,  +  (w,  (c^A,-,  +  fc^A.) 
4A,A,_i  -  ru>,(A, -r  A,_i) 


with  c\0)  =■-  c,.  The  user  can  input  the  number  of  iterations  and  the  convergence  factor  (.  These 
two  parameters  allow  some  flexibility  in  how  much  debiasing  Ray  applies. 

On  each  iteration  the  sound  speed  at  any  depth  is  only  affected  by  its  two  nearest  neighbors. 
Thus  the  number  of  iterations  controls  how  local  or  global  the  debiasing  will  be.  The  convergence 


factor  controls  how  close  Ray  tries  to  get  to  the  original  sound  speeds.  If  it  is  set  to  1.0,  then  Ray 
will  attempt  to  compensate  for  all  of  the  smoothing.  If  it  is  set  to  0.5  then  Ray  will  only  try  to  get 
rid  of  one  half  of  the  offset  caused  by  the  smoothing.  For  a  few  test  profiles,  ten  iterations  with  e 
set  to  1.0  converges  to  the  original  sound  speeds. 

An  example  of  a  smoothed  sound-speed  profile  with  the  bias  removed  is  displayed  in  Figure  1. 


Figure  1.  The  solid  line  shows  a  section  of  a  piecewise  parabolic  sound  speeu  profile  generated  by  Ray.  The 
circles  show  the  input  sound  speeds  which  were  used  by  Ray  to  generate  the  profile.  The  curve  to  the  right 
shows  the  piecewise  constant  second  derivative  of  the  solid-line  profile,  with  sections  centered  on  the  input 
depths.  In  general  there  can  be  linear  sections  between  the  parabolic  sections  but  none  are  shown  in  this 
example.  The  dashed  linear  segments  show  the  constant  gradient  approximation  used  in  circular  arc  ray 
tracing  programs  such  as  MPP  and  RDRYT. 
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Range  Dependence 

The  way  that  Ray  handles  range  dependence  reflects  the  design  philosophy  of  creating  a  simple 
model,  consistent  with  the  input  data,  that  produces  smooth  wavefronts.  Sound  speeds  are  input 
as  tables  of  sound  speed  versus  depth  at  increasing  ranges  from  the  source.  At  each  range  the  sound 
speed  is  smoothed  and  debiased  as  described  above.  Three  different  models  of  range  dependence 
are  provided. 

The  simplest  model  treats  the  ocean  as  a  series  of  range  independent  sections.  At  each  range 
where  a  sound-speed  profile  is  defined  the  entire  profile  jumps  abruptly  to  the  new  values.  No 
compensation  is  made  for  obeying  Snell’s  law  at  the  interface.  The  rays  that  make  up  a  wavefront 
pass  through  each  interface  at  different  depths.  Since  the  amount  the  sound  speed  jumps  varies 
with  depth,  this  model  introduces  a  modulation  of  the  wavefront  for  each  interface  that  is  passed 
through.  There  are  some  regions  of  the  ocean  where  this  modulation  does  not  produce  serious 
problems  but  there  are  many  regions  where  it  is  unacceptable.  All  attempts  to  eliminate  this 
modulation  by  making  some  correction  at  the  interface,  such  as  obeying  Snell’s  law,  have  failed. 
The  corrections  have  changed  the  modulation  but  have  not  eliminated  it. 

The  next  level  of  sophistication  lets  the  sound-speed  profile  vary  linearly  with  range.  As¬ 
sume  two  adjacent  profiles  have  been  input  at  ranges  rj  and  rJ+i,  then  the  sound  speed  at  any 
intermediate  range  r  will  be 

c(?,r)  =  r(r;,z)+  (c(rJ+l,z)-c(rJ,z)). 

r;+t  “  rj 

The  equations  of  motion  (see  above)  for  a  ray  depend  upon  c,  dc/dz  and  dc/dr.  In  this  model  c  and 
dc/dz  vary  smoothly  with  range  and  depth  while  dc/dr  changes  discontinuously  at  each  interface. 
If  this  term  is  important  in  determining  the  path  of  a  ray  then  the  linear  range  dependent  model 
will  also  contain  modulations  of  the  wavefront  that  depend  upon  the  depths  at  which  each  ray 
passes  through  the  interfaces. 

Ray  prov’das  two  implementations  of  the  linear  range  dependent  model  in  order  to  establish 
whether  the  contributions  due  to  the  dc/dr  term  are  significant.  One  of  them  keep6  this  term 
and  the  other  one  drops  it.  Using  input  from  the  Levitus  data  base,  the  differences  between  the 
wavefronts  generated  by  these  two  implementations  have  been  of  the  same  scale  as  the  numerical 
noise.  This  evidence  has  led  us  to  conclude  that  the  linear  range  dependent  model  is  sufficient  for 
the  intended  use  of  Ray. 

If  Ray  is  run  with  an  environment  where  dropping  the  dc/dr  term  produces  a  significant  change 
in  the  wavefront  then  it  is  very  likely  that  the  linear  range  dependent  model  is  inadequate  for  that 
environment  even  if  the  dc/dr  term  is  included. 
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Bathymetry  and  Bouncing 

Version  1.0  of  Ray  provides  several  different  options  for  dealing  with  bathymetry.  The  soft 
option  keeps  a  record  of  the  smallest  distance  to  the  bottom  for  each  ray  of  a  wavefront.  Rays 
are  still  traced  even  if  they  go  below  the  bottom.  In  these  cases  the  distance-to-bottom  becomes 
negative  and  the  most  negative  value  is  recorded.  The  range  where  the  smallest  (or  most  negative) 
distance-to-bottom  occurs  is  also  recorded.  The  absorbing  option  terminates  any  ray  that  touches 
the  bottom.  The  reflecting  option  reflects  rays  that  hit  the  bottom. 

The  bathymetry  data  that  Ray  reads  in  can  be  smoothed  in  the  same  manner  that  the  sound 
speeds  are  smoothed.  Only  one  width  parameter,  bathjaoothing,  may  be  input.  The  delault  value 
is  10  km.  It  is  used  for  all  bathymetry  points  whose  nearest  neighbor  is  at  least  4  times  greater 
than  this  width  away.  All  bathymetry  points  that  have  a  nearest  neighbor  closer  than  4  times  the 
bathymetry  smoothing  width  have  their  smoothing  width  set  to  one  quarter  of  the  distance  to  the 
nearest  neighbor.  This  scheme  creates  alternating  sections  of  straight  lines  and  parabolas.  If  the 
bathymetry  smoothing  width  is  set  to  0  then  no  smoothing  of  the  bathymetry  is  performed  and 
Ray  uses  a  bottom  made  up  of  straight  lines  connecting  the  input  bathymetry  points.  No  attempt 
has  been  made  to  “debias”  the  bathymetry,  although  the  bathymetry  points  can  be  debiased  before 
being  passed  on  to  Ray. 

The  heart  of  a  bouncing  algorithm  consists  of  finding  exactly  where  a  ray  path  first  intersects  a 
boundary.  Ray  uses  the  following  approach  to  find  these  intersection  points.  For  every  integration 
step  taken,  a  quick  check  is  performed  to  determine  if  the  ray  might  have  possibly  crossed  a 
boundary.  Let  n  and  r3  denote  the  range  at  the  beginning  and  end  of  a  step.  If  the  possibility  of 
crossing  exists  then  the  depth  of  this  short  section  of  the  ray  path  and  the  depth  of  the  boundary 
are  both  parameterized  as  parabolic  in  range  and  the  difference  of  these  two  parabolas  is  taken  to 
give  the  distance  from  the  bottom  as  a  function  dr  =  rj  -  r, 

^r»y  —  ^bottom  ~  ^mi» 

as  zq  +  dr  z\  +  dr2  z3. 

We  assume  that  if  the  integration  step  size  has  converged  then  the  first  crossing  of  the  boundary 
will  occur  either  near  the  smallest  positive  real  root  of  this  equation  or  near  the  minimum  of 
this  equation  if  it  has  no  real  roots.  In  either  case  we  take  one  integration  step  to  the  range 
rf1)  =  ri  +  r^i,t  which  is  a  rough  guess  of  where  the  first  intersection  will  occur.  The  expansion 
above  is  then  repeated,  this  time  expanding  about  r^.  If  no  real  roots  are  found  here  then  we 
assume  that  the  ray  does  not  cross  the  boundary  at  this  step.  If  real  roots  are  found  then  our 
next  guess  for  the  range  where  the  first  intersection  occurs  is  =  r^  +  r^t,  which  is  found  by 
adding  the  root  of  the  new  parameterization  with  the  smallest  absolute  value  to  r^.  This  process 
is  repeated  until  either  no  real  roots  are  found  between  ri  and  r3  or  <  z’tol-  The  default 
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value  for  rtoj  is  10_s  m.  If  no  roots  are  found  then  we  did  not  hit  the  boundary  in  this  step.  If  the 
tolerance  condition  is  met  in  the  ith  iteration  then  we  consider  that  the  ray  hit  the  boundary  in 
this  step  at  the  range  and  the  appropriate  action  is  taken  depending  on  the  type  of  bathymetry 
requested . 
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3.  Program  Input 

Invoking  Ray  with  no  arguments  will  produce  an  output  similar  to 

/ft************************  Ray  ♦*•***••**•****•••••••••*/ 

/• 

•  Version  1.00  Friday,  loveaber  13,  1902,  3:33  pa 

a  Copyright  (C)  1992,  Woods  Hols  Oceanographic  Institution. 

•  Use  -license  flag  for  use,  copying  and  distribution  conditions. 

•/ 

usage: 


ray 

initfile.ray 

[cBtndl]  Ccand2]  . . .  [eandl]  (flags] 

•here : 

initfile.ray 

is  a  ray  Initialization  fils 

caadl 

. . .  candf 

are  cosnand  line  paraaeters  used  in  the 

initialization  file  via  tit  ...  tit 

flags: 

-d 

(debug) 

send  sons  debugging  info  to  stdout 

-1 

(license) 

display  the  license  agrseaent 

-■ 

(aakeinit) 

generate  an  initialization  fils 

-h 

(help) 

send  detailed  help  to  stdout 

~P 

(parse) 

send  parsed  lnit  fils  to  stdout 

-¥ 

(verbose) 

send  information  to  stdout  as  ve  work 

which  lists  all  of  Ray’s  command  line  options.  Running  Ray  with  the  -aakeinlt  option  as  in 

Ray  -a 

causes  Ray  to  print  a  prototype  initialization  file  on  the  standard  output.  Capturing  and  editing 
this  file  is  an  easy  way  to  quickly  operate  Ray. 

The  initialization  file  can  then  be  sent  back  into  Ray  by  specifying  its  name  as  the  first  command 
line  parameter  as  in 

Ray  initialize. ray  , 

which  will  run  the  Ray  program  with  “initialize. ray”  as  the  initialization  file.  It  is  always  the  first 
file  read  by  Ray.  Its  format,  structure,  and  syntax  are  described  in  this  section.  The  formats  of  the 
other  input  files  are  described  in  an  appendix. 
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Initialization  File  Format 


The  initialization  file  format  will  be  familiar  to  people  who  work  with  the  C  programming 
language.  All  text  enclosed  between  /•...*/  is  a  comment  and  is  ignored.  All  of  the  initialization 
parameters  are  organized  into  groups.  The  parameters  within  a  group  can  be  specified  in  two 
different  ways.  One  way  is  a  single  line  (or  statement)  for  each  parameter.  For  example, 
model  integration  =  rk.2; 
model  range.depend  ■  grad-z; 

will  set  the  integration  routine  to  be  and  the  range  dependence  to  be  of  the  type  grad_z.  More 
information  on  exactly  what  these  mean  will  follow  below.  Notice  the  semicolon  at  the  end  of  each 
statement.  The  second  method  of  specifying  parameters  within  a  group, 
model  { 

integration  =  rk.2; 
rangm-dmpmnd  =  grad-2; 

>; 

sets  these  two  parameters  exactly  like  the  example  given  previously.  This  method  saves  on  redun¬ 
dant  typing  and  encourages  users  to  put  all  of  the  parameters  within  a  single  group  together.  There 
are  currently  six  types  of  parameters  that  may  be  specified:  string,  dimensioned  numerical,  array 
of  dimensioned  numerical,  dimensionless  numerical,  choice,  and  flag. 

String  Parameters 

String  parameters  are  used  for  specifying  filenames.  The  name  must  always  be  inside  a  pair  of 
double  quotes,  as  in 

input  prof. film  »  "tmmtl.mmp"; 

Dimensioned  Numerical  Parameters 

Dimensioned  numerical  parameters  are  input  as  a  number  followed  by  units  in  parentheses. 
For  example 

rmcmivmr  range  ■  1000.0  (km); 

specifies  that  the  receiver  is  at  a  range  of  1,000  kilometers.  The  inclusion  of  units  is  not  optional 
but  choices  of  units  are  available.  Two  kinds  of  units  are  used,  lengths  and  angles.  Internally,  all 
angles  are  in  radians  and  all  lengths  are  in  meters.  Table  1  lists  valid  length  and  angular  units. 
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Table  1.  Valid  length  and  angular  units. 


unit 

name 

conversion  to  meters 

(n) 

meter 

l.OeO 

<kn) 

kilometer 

1.0e3 

(Ini) 

nautical  mile 

1.852e3 

(ft) 

foot 

3.04 8e-l 

(furlong) 

furlong 

2.01 17e2 

(parsec) 

parsec 

3.804el6 

(cm) 

centimeter 

1.0e-2 

(nn) 

millimeter 

1.0e-3 

(1(a) 

megameter 

1.0e6 

(lightyear) 

lightyear 

9.46el5 

(angstroa) 

angstrom 

1.0e-10 

unit 

name 

conversion  to  radians 

(radians) 

radians 

l.OeO 

(degrees) 

degrees 

3  14159265358979  /  180.0 

Array  of  Dimensioned  Numerical  Parameters 

An  array  of  dimensioned  numerical  parameters  is  specified  like  a  single  dimensioned  numerical 
parameter  except  the  single  number  is  replaced  by  a  left  brace,  a  list  of  numbers  and  a  closing  right 
brace  as  in 

angles  { 

specific  « 

{1.0,  2.0,  3.0,  4.0,  S.O  >  (degrees); 

> 

Note  that  the  list  of  numbers  can  be  delimited  by  either  whitespace  characters  (which  consist  of 
space,  tab,  linefeed,  and  carriage  return)  or  commas  or  both. 

Dimensionless  Numerical  Parameters 

Dimensionless  numerical  parameters  are  similar  to  dimensioned  numerical  parameters  but  they 
must  not  have  any  unit  specification,  even  an  empty  “()"  is  not  allowed.  For  example, 
nodal  Margins  »  100; 
is  a  valid  statement. 


Choice  Parameters 

Choice  parameters  give  the  user  a  specific  range  of  choices  such  as 
nodal  integration  ■  rk_2; 

which  will  set  the  integration  routine  to  be  rL2.  Other  choices  for  this  parameter  are  rk_23  and 
rk.4.  If  a  choice  parameter  is  misspelled  or  otherwise  inappropriate,  then  an  error  message  will  be 
generated  showing  where  the  problem  occurred  and  what  the  valid  choices  are. 
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Flag  Parameters 

Flag  parameters  look  similar  to  choice  parameters.  A  flag  parameter  is  set  by  mentioning  its 
name  as  in: 

output  uavsfront; 

which  will  ensure  that  all  of  the  wavefront  variables  for  this  run  are  put  into  the  output  file. 

Initialization  Parameters  by  Group 

There  are  seven  groups  of  parameters  that  are  specified  in  the  initialization  file.  They  are: 
input,  output,  .ourco,  recti vor,  angles,  paths,  and  aodel.  All  of  the  initialization  parameters  are 
described  below,  ordered  by  group. 

Input  Group 

The  input  group  specifies  the  additional  files  that  will  be  read  in  by  Ray.  For  example 
input  { 

protJile  =  "teatl.ssp"; 
bath Jilt  •  "teatl.btb"; 

>; 

specifies  that  the  file  “testl.ssp”  will  be  read  in  as  a  profile  file  containing  sound-speed  profiles,  and 
the  file  “testl.bth”  will  be  read  in  as  a  bathymetry  file  containing  bathymetry.  In  order  to  maintain 
compatibility  with  existing  software,  the  profiles  and  bathymetry  may  be  specified  together  in  a 
single  MPP  file  with  the  line 

input  app.fila  *  "taatl.npp"; 

If  an  nppjrilt  is  specified,  then  it  is  an  error  to  specify  either  a  prof_fil*  or  a  bath_f  il«.  The 
profile  and  bathymetry  file  formats  are  detailed  in  appendices.  The  npp.1  ila  option  is  only  included 
to  maintain  compatibility  with  existing  software.  This  file  format  is  not  recommended  and  is  not 
described  in  this  report. 

Output  Group 

The  output  group  specifies  the  name  of  the  output  file  generated  by  Ray,  and  what  variables 
will  be  included  in  this  file.  The  line 

output  aat.file  ■  "tMtl.aat"; 

tells  Ray  to  name  the  output  file  “testl.mat”.  The  user  can  tailor  what  will  be  included  in  the 
output  file  with  the  following  flag  parameters, 
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output 

effect 

initialization; 

save  initialization  parameters 

tilanaaas; 

save  filenames 

data; 

save  the  time  and  date  of  the  run 

»ound_apaad« ; 

save  the  sound-speed  tables 

bathyaatry; 

save  the  bathymetry  table 

paths ; 

save  path  information  along  each  ray 

■aval ront ; 

save  all  ray  parameters  at  receiver  range 

■vary thing; 

save  all  of  the  above 

anvlrouant-only ; 

don’t  trace  rays 

If  no  output  flags  are  specified  then  Ray  will  assume  everything  should  be  saved.  To  save  the  paths 
the  user  must  explicitly  set  a  paths  fixod^dr  or  a  paths  stapa_p«r  in  addition  to  the  output  flag. 
This  is  a  safety  feature  to  prevent  accidentally  filling  a  disk  with  detailed  path  information  from 
thousands  of  rays.  If  onvironnont-only  is  chosen  then  no  rays  will  be  traced  regardless  of  the  state 
of  the  other  output  flags. 

5ounce  Group 

The  source  group  specifies  the  depth  and  range  of  the  source.  The  lines 
source  { 

dspth  «  600.00  (■); 
rang*  ■  1000  On) ; 

> 

will  put  the  source  600  m  below  the  surface  at  a  range  of  1,000  km 
Receiver  Group 

The  receiver  group  specifies  the  depth  and  range  of  the  receiver.  The  depth  of  the  receiver  is 
not  used  by  Ray  Version  1.0.  For  example, 
r«c«iv«r  < 

depth  ■  1000.0  (■); 
rang*  •  3000 . 0  (ka) ; 

> 

specifies  a  receiver  at  a  depth  of  1000  m  and  a  range  of  3000  km. 

Angle  Group 

The  angle  group  contains  all  of  the  information  needed  to  specify  the  launch  angles  of  all  the 
rays  that  will  be  traced.  There  are  two  ways  to  specify  these  angles,  one  is  to  provide  two  extreme 
angles  and  the  total  number  of  (equal  spaced)  angles  as  in 


19 


angles  { 

first  *  IS  (dsgrsss); 
last  *  -16  (dsgrsas); 
nuabar  *  3; 

> 


which  tells  Ray  to  6hoot  rays  at  the  angles  15,  0,  and  -15  degrees.  The  user  may  alternatively 
specify  specific  angles  to  shoot  as  in 
angles  specific  »  { 

1.0, 

2.0. 

3.0, 

> (degrees)  ; 


which  tells  Ray  to  shoot  three  rays  at  initial  angles  of  1,  2,  and  3  degrees.  It  is  an  error  to  specify 
angles  in  both  formats. 


Paths  Group 


The  parameters  in  the  paths  group  allow  the  user  to  save  information  about  each  ray  as  it  is 
traced.  For  example 
paths  { 

■izurange  >0.0  (taa) 
aajLxange  *  1000  On) 
fixed-dr  ■  i  (ka) 


will  save  path  information  for  every  ray  as  it  is  traced,  every  km  over  the  region  from  0  to  1,000 
km.  If  the  aimrange  and  aax-range  are  not  specified  (or  if  they  are  set  to  zero)  then  the  path 
region  is  set  to  the  entire  region  between  the  source  and  the  receiver.  If  stsps.p*r  is  used  instead 
of  f  ixa<Ldr  then  Ray  will  output  information  along  each  ray  every  staps.par  steps.  Since  the  step 
size  can  vary  with  depth,  the  range  spacing  of  each  path  will  vary.  The  information  that  is  stored 
at  points  along  a  path  can  be  tailored  by  specifying  the  following  paths  colons  flags: 
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paths  colunns 

contents 

range; 

range  of  the  point 

depth ; 

depth  of  the  point 

tine; 

time  it  took  to  get  to  the  point 

angle; 

angle  of  the  ray  at  the  point 

speed ; 

speed  of  sound  at  the  point 

grad; 

d  c/8z  at  the  point 

top.bouaeea ; 

number  of  top  bounces 

bo  t  .bounce* ; 

number  of  bottom  bounces 

ray  jmnbar ; 

which  ray  the  point  belongs  to 

every thing; 

all  of  the  above 

If  no  columns  are  specified  then  all  the  columns  are  used.  In  addition  if 
paths  include  bounces; 

is  specified  then  every  point  where  a  ray  bounces  is  included  in  the  path  in  addition  to  all  of  the 
points  generated  by  the  iixad_dr  specification  or  the  stepa.psr  specification.  In  order  to  only  save 
the  bounce  points  set  iixad-dr  to  a  range  that  is  greater  than  the  separation  between  the  source 
and  the  receiver. 

Model  Group 

The  model  group  of  parameters  specify  how  Ray  will  model  the  ocean.  The  nodal  integration 
choice  controls  which  integration  routine  is  used.  The  options  are 


nodal  integration  ■  effect 


rkJt;  2nd  order  Runge  Kutta 

rk-23;  2nd  -  3rd  order  Runge  Kutta 

rk.4;  4th  order  Runge  Kutta 

The  nodal  raagn.dspand  choice  controls  how  Ray  interpolates  between  different  sound-speed  pro¬ 
files.  There  sue  three  options;  . 


nodal  rangejdapead  *  effect 


none; 

gradj; 

full; 


no  interpolation 

interpolate  c  and  dc/dz  but  ignore  dc/dr 
interpolate  c  and  dc/dz  and  include  dc/dr 


The  nodal  bathynatry  ■  choice  controls  how  Ray  deals  with  bathymetry. 
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■odel  bathymetry 


effect 


non*;  flat  bottom,  no  bathymetry 

•oft;  pass  through  bottom 

absorbing;  stop  at  bottom 

reflecting;  reflect  off  bottom 

The  soft  choice  records  the  closest  approach  to  the  bottom  (or  the  deepest  descent  into  the  bottom) 
for  each  ray.  It  also  records  the  range  where  this  closest  approach  occurs.  The  absorbing  choice 
terminates  any  ray  that  touches  the  bottom.  The  reflecting  choice  lets  rays  bounce  off  of  the 
bottom. 

Within  the  aodel  group  is  a  subgroup  that  controls  the  size  of  the  range  step  taken  by  the 
integration  routines.  Ray  provides  a  facility  for  a  depth  dependent  step  size.  The  shape  of  this 
function  is  determined  by  the  prof-saoothlng  widths  which  are  described  below.  The  overall  size 
of  the  steps  can  be  controlled  by  setting  the  step  size  multiplier,  such  as 
■odel  rang* _atep  nultiplier  *  O.S; 

which  halves  all  of  the  range  steps.  The  maximum  and  minimum  range  step  can  also  be  controlled. 
This  process  occurs  after  the  multiplier  has  been  applied  so  that  the  units  used  for  the  minimum 
and  maximum  are  not  scaled  by  the  multiplier.  The  lines 
■odel  range-»tep  { 
a ax  >  10  (a); 

■in  ■  10  (■); 

> 

will  set  the  range  step  size  to  10  m  for  all  depths  regardless  of  the  value  of  the  multiplier.  The 
maximum  and  minimum  can  also  be  set  to  different  values.  If  the  maximum  is  less  than  the 
minimum,  then  the  step  size  is  set  to  the  minimum. 

There  are  two  parameters  in  the  aodel  group  that  control  how  dose  Ray  stays  to  the  input 
sound-speed  profiles.  These  are  dimensionless  numerical  parameters  and  are  set  as  in  the  following 
example 

■odel  deblaa  { 
factor  ■  1.0; 
iteration  ■  10; 

> 

A  detailed  description  of  what  these  two  parameters  do  is  contained  in  the  Computational  Algorithm 
section  of  this  report. 

There  are  six  numerical  parameters  in  the  aodel  group,  these  parameters  and  their  default 
values  are: 

■odel  { 

prof  jaoothing  ■  { 

6,  6,  S,  5,  S,  20,  12,  13,  12,  26,  26.  26.  60, 
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60,  60,  60,  60,  60,  60,  60,  60,  60,  60,  60.  60, 

126,  126,  260,  260,  260,  260,  260,  260,  260 
>(*); 

bath.saoothi.Ag  *  10  (ka) ; 
bottom-depth  *  6000  (a); 
earth-radius  *  6378.137  (a); 

tolerance  »  la-06  (a); 
aarglB*  *  100  ; 

> 

The  prof-smoothing  array  contains  the  smoothing  widths  for  the  sound-speed  profiles  as  discussed  in 
the  Algorithms  section.  The  bath.aaoothing  parameter  contains  the  single  width  used  for  smoothing 
the  bathymetry.  The  bottomjdepth  parameter  controls  the  absolute  maximum  depth  that  can 
be  used  in  a  ray  trace.  If  there  is  bathymetry  specified  that  is  deeper  than  bottoa.depth  then 
bottoajdspth  is  automatically  set  to  the  deepest  bathymetry  point.  The  earth-radius  parameter  is 
used  for  the  spherical  earth  correction.  If  it  is  set  to  0  then  no  spherical  earth  correction  is  used. 
The  z.tolerance  parameter  controls  how  close  a  ray  need  get  to  a  boundary  before  it  bounces  off 
of  it  as  described  in  the  Bathymetry  and  Bouncing  section.  It  is  currently  set  to  10“6  m.  The 
margin*  parameter  controls  how  much  extra  depth  is  allowed  over  the  top  and  under  the  bottom. 

The  default  values  of  these  numerical  parameters  should  work  fine  in  most  situations.  If  they 
are  set  incorrectly  spurious  results  may  be  obtained. 

Command  Line  Substitution 

Often  the  situation  arises  when  several  ray  traces  need  to  be  performed  which  differ  in  only 
one  or  two  parameters.  To  facilitate  work  on  such  problems,  Ray  has  the  ability  to  do  command 
line  substitution  in  the  initialization  file.  This  allows  repeated  use  of  the  same  initialization  file 
for  a  number  of  different  ray  traces.  Everywhere  in  the  initialization  file  where  a  $1$  occurs,  the 
second  command  line  parameter  is  substituted  for  the  tit.  Everywhere  that  a  t2t  is  encountered, 
the  third  command  line  parameter  '13  substituted  for  every  t2t,  and  so  on.  For  example  suppose 
the  “init.ray"  contains  the  lines 

input  prof-file  *  "tit. sap"; 
input  bath-file  *  "lll.bth"; 
nodal  integration  ■  t2t; 
and  then  suppose  that  Ray  is  invoked  by 
Ray  init .ray  run 17  grad-Z 

then  this  will  have  the  same  effect  as  if  the  initialization  file  contained 
input  prof  .file  ■  "runl7.sap”; 
input  bath-file  «  "runl7.bth"; 

■odel  integration  ■  gred_* ; 

Note  that  although  replaceable  parameters  are  allowed  to  split  a  string  (as  in  "$l$.prf"  above) 
it  is  not  valid  to  try  to  split  up  a  word  that  is  not  in  quotes  with  a  replaceable  parameter.  For 
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example,  the  liue 

model  integration  rk.Slt; 
is  not  valid  and  will  generate  several  error  messages. 

Although  Ray  has  a  large  number  of  inputs  parameters  which  make  it  adaptable  to  a  variety  of 
needs,  most  users  will  never  want  to  access  all  of  these  parameters  since  Ray  has  a  reasonable  set  of 
default  values.  However,  There  is  some  information  that  needs  to  be  included  in  the  initialization 
file  in  order  for  Ray  to  run.  The  following  listing  shows  a  very  short  initialization  file  which  contains 
all  of  the  needed  information. 

/•*••  Minimum  information  initialization  file  ***/ 


input  prof .file 
output  mat.file 
•ourca  dopth 
racaivar  range 
anglaa  firat 
anglas  last 
anglaa  number 


teetl.eap"; 

testl.mat"; 

600  (m); 

220  (km); 

IS  (degrees) ; 
-15  (degrees); 
500; 


If  a  bathymetric  ray  trace  is  desired  then  add  the  line 
input  bath.fila  *  "fnams.bth"; 

and  Ray  will  read  in  the  bathymetry  information  from  this  file  and  set  modal  bathymetry  *  re¬ 
ft  acting; .  If  path  information  is  desired  then  add  a  line  such  as 
paths  fixsd_dr  =  1000  (m); 

which  will  cause  Ray  to  save  information  along  the  path  of  each  ray  every  1000  m. 
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4.  Program  Operation 

Ray  is  written  in  the  ANSI  C  language.  The  source  code  io  split  up  into  several  files  (which 
are  called  "translation  units’"  in  C  jargon).  All  of  these  fiies  have  a  “.c’’  extension.  Many  have  an 
associated  header  file  with  the  same  name  and  a  u.h”  extension.  A  general  outline  of  what  Ray  does 
when  it  runs  follows  below.  Each  major  section  of  the  outline  includes  the  name  of  the  translation 
unit  which  is  primarily  in  charge  of  that  section. 

Outline 

I.  Initialization:  initray.e 

a.  Print  banner  and  copywrite  notice. 

b.  Read  command  line  parameters.  Exit  here  if  there  are  unknown  command  line  parameters. 

c.  Open  and  read  initialization  file  (there  are  errors  in  initialization  file  then  print  error 
messages  and  exit). 

d.  Open  appropriate  MPP,  profile  and  bathymetry  files.  Open  the  output  file. 

e.  Save  all  initialization  parameters,  start  time,  filenames  in  the  output  file. 

II.  Read  auxiliary  input  files:  mppio.c 

a.  Read  MPP  file  containing  sound-speed  profiles  and  bathymetry. 

b.  Read  profiles  file. 

c.  Read  bathymetry  file. 

III.  Pre-process  the  environment:  preray.c 

a.  Get  maximum  number  of  sound-  speed  depths. 

b.  Make  a  table  of  widths. 

c.  Make  idex  tab1  . 

d.  Make  rstep  table. 

e.  Debias  the  sound  speeds. 

f.  Make  sound-speed  tables. 

g.  Make  speed  at  receiver  table. 

h.  Make  bathymetry  table. 

j.  Save  environment  (if  requested  to  do  so). 

IV'.  Do  the  ray  trace:  ray.c 

a.  For  every  angle:  shoot  one  ray. 

b.  Save  wavefront  information. 

Program  Function  by  Translation  Unit 

helpray.c  Contains  the  text  of  Ray’s  extensive  help  documentation  and  a  little  function  to  print  it 
out. 

initray.e  Prints  out  the  banner  and  copywrite  m  ;ce.  Reads  ana  per  cesses  command  line  parame¬ 
ters.  If  any  errors  occur  (such  as  an  unrecognized  command  line  option)  the  usage  text  is 
displayed  along  with  an  error  message.  Reads  and  processes  the  initialization  file.  If  any 
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errors  occur  while  reading  the  initialization  hie,  each  offending  line  is  printed  out  along 
with  an  error  message  showing  where  in  the  line  the  error  occurred  and  what  the  error  was. 
The  total  number  of  errors  is  displayed  and  then  the  program  exits.  If  the  initialization 
hie  is  parsed  successfully  then  its  contents  are  checked  for  ambiguities.  If  any  ambiguities 
occurred  then  all  ambiguities  in  the  initialization  file  are  reported  and  the  program  exits. 
If  no  ambiguities  occurred  then  we  open  the  other  input  files  and  the  output  file, 
mppio.c  Reads  MPP,  profile  and  bathymetry  files.  Puts  all  of  this  information  into  a  single  app 
structure.  Since  we  use  reada.c  (see  below)  to  read  ASCII  files,  comments  may  be  freely 
placed  throughout  these  files  be  inclosing  them  in  /•  . .  .*/. 
outray.c  Contains  functions  to  write  variables  to  the  output  file  in  binary  MatLab  format.  Ray  has 
one  machine  dependent  parameter.  It  is  called  Matjaach-nua  and  is  located  in  outray.h. 
There  is  a  nearby  comment  which  explains  what  to  set  this  parameter  to  according  to 
which  machine  Ray  will  be  compiled  on.  If  it  is  set  incorrectly,  Ray  will  run  without  errors 
but  MatLab  will  be  unable  to  read  Ray’s  output. 

preray.c  Preprocesses  the  environment.  Constructs  an  index  table  and  tables  of  sound-speed  param- 
eterizations,  bathymetry  parameterizations  and  automatic  rangestep  parameterizations. 
All  of  these  functions  are  accessed  through  the  function  preprocessi). 
ray.c  Does  the  actual  1  ng.  Contains  all  of  the  functions  for  integrating  the  equations 

of  motion  and  bounc.ng  rays  off  of  the  surface  and  the  (flat)  bottom.  Also  contains  the 
function  mcin()  which  is  the  main  routine  for  ray. 
reada.c  Reads  ascii  files  and  splits  them  into  tokens.  A  token  is  defined  as  a  contiguous  series  of  non¬ 
whitespace  characters  separated  by  whitespace  characters  (space,  tab,  linefeed,  and  carriage 
return)  and  delimiter  characters  (  *,/"  {}();#).  Each  delimiter  character  is  a 
token.  Strings  and  comments  are  handled  appropriately.  All  text  between  a  /•  and  a 
*/  is  ignored.  All  text  between  pairs  of  double  quotes  (  "  . . .  ")  is  treated  as  a  string. 

Comments  may  be  nested. 

utils.c  Contains  small  functions  and  macros  used  by  many  of  the  other  translation  units.  The 
macros  are  in  the  file  utils.h. 
version.c  Contains  the  version  number  and  version  date. 

Detailed  Outline  of  Tracing  One  Ray 

Ray  uses  two  different  structures  for  tracing  rays.  The  structure  rlnfos  contains  all  of  the 
information  about  a  ray  that  is  used  in  creating  a  wavefront.  The  icllowing  fragment  from  ray.h 
details  all  of  the  components  of  the  rlnfos  structure. 


struct 

rlnfos  < 

/* 

all  start,  end  info  about  a  single  ray 

*/ 

double 

s angle, 

/* 

angle  of  ray  at  source 

*/ 

sdsptb, 

/* 

depth  of  ray  at  source 

*/ 

srangs, 

/* 

range  of  the  source 

*/ 

r angle, 

/* 

angle  of  ray  at  reviver 

*/ 

rdepth. 

/* 

depth  of  ray  at  rscaivsr 

♦/ 

rrangs. 

/* 

range  of  receiver 

a/ 

tins. 

/* 

tine  of  arrival  at  the  receiver 

*/ 

top.turas , 

/- 

nuaber  of  top  tarns 

•/ 

bot.turns , 

/* 

nanber  of  botton  tarns 

•/ 

top.bonks , 

/* 

nuaber  of  top  bounces 

*/ 

bot.bonks , 

/* 

nanber  of  botton  bounces 

•/ 

tot.turns , 

/* 

total  nanber  of  tarns 

*/ 

signat, 

/• 

■  10000  •  bonks  ♦  turns 

•/ 

utp. 

/• 

upper  turning  point  at  receiver 

*/ 

ltp. 

/• 

lover  turning  point  at  receiver 

*/ 

aiss. 

/• 

niss  flag 

•/ 

arrival. 

/• 

arrival  type 

*/ 

c_at_s. 

/• 

speed  of  sound  at  source 

*/ 

nadir, 

/• 

closest  approach  to  botton 

•/ 

nadir. at ; 

/• 

range  of  closest  approach  to  botton 

•/ 

rinfos 

•prsv; 

/• 

not  used 

•/ 

r infos 

•next ; 

/• 

points  to  next  rinfo 

a/ 

>; 

The  other  structure  used  for  tracing  rays  is  much  smaller.  It  is  called  rels  which  stands  for  ray 
elements.  The  components  of  this  structure  are  what  actually  get  integrated  when  a  ray  is  being 
traced.  The  components  of  rels  and  their  significance  are  given  in  following  code  fragment  which 
is  also  from  ray.h 


struct  rels  { 

/*  this  is  the  structure  vs  integrate  •/ 

double  r, 

/*  range  a/ 

z. 

/•  depth  a/ 

s. 

/•  sin(thata)  a/ 

e. 

/•  cos(theta)  */ 

t, 

/•  tine  •/ 

a; 

/*  theta:  only  used  in  derivitives,  *0T  integrated  a/ 

>; 

Note  that  0,  the  angle  of  the  ray  is  not  integrated.  Instead,  cos0  and  sin  9  are  integrated  separately. 
This  allows  the  integration  to  be  performed  with  only  addition,  subtraction,  multiplication,  and 
division  operations.  No  calls  to  transcendental  functions  are  needed. 

As  described  in  the  Algorithms  section,  the  slice  of  the  ocean  we  are  tracing  rays  through  is 
broken  up  into  a  sequence  of  regions  such  that  the  sound  speed  in  each  region  varys  linearly  with 
range.  If  the  source  is  located  before  the  first  profile,  or  if  the  receiver  is  located  after  the  last 
profile  then  profiles  are  repeated  as  necessary  in  order  to  assure  that  the  entire  path  of  a  ray  is 
bounded  between  profiles. 
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For  each  region,  global  pointers  to  the  proper  sound  speed  tables  sure  set  and  then  the  ray  is 
shot  through  the  region  by  repeating  the  following  steps: 

1)  compute  the  range  stepsize 

2)  take  one  Runge  Kutta  step 

3)  do  quick  checks  of  boundary  crossing. 

4)  If  a  boundary  might  have  been  crossed  then  do  exact  test  for  a  boundary  crossing. 

5)  If  a  boundary  was  crossed  take  the  appropriate  action  depending  on  the  setting  of  aodel 
bathymetry, 

6)  Count  turning  points. 

When  the  ray  would  get  to  within  one  meter  of  the  end  of  a  region  the  range  stepsize  is  set  so  that 
the  ray  will  land  right  on  the  boundary  of  the  region.  If  there  is  another  region  to  go  through,  then 
we  are  ready  to  repeat  the  above  process,  otherwise  we  are  at  the  receiver  range  and  the  tracing  of 
this  ray  is  finished. 
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5.  Program  Output 

The  output  from  Ray  is  a  single  Matlab  file  ( Moler ,  Little  and  Bangert ,  1987].  The  name  of 
this  file  is  specified  in  the  initialization  file  on  a  line  of  the  form 
output  mat-fils  *  "foms  .mat" ; . 

This  file  is  written  in  binary  double  precision  MatLab  format  and  can  be  read  directly  by  the 
MatLab  program  on  a  variety  of  platforms.  It  contains  a  collection  of  MatLab  variables  that  can 
totally  specify  the  inputs  and  outputs  of  the  Ray  program.  The  variables  are  divided  up  into  seven 
groups.  The  user  has  control  over  which  group  or  groups  of  variables  get  saved  by  specifying  output 
groupjsaae;  in  the  initialization  file.  The  statement  output  everything  causes  Ray  to  do  just  that. 

The  names  of  the  seven  output  groups  are  given  in  the  table  below  along  with  a  brief  description 
of  what  each  group  contains. 

output  group  contents 


initial izatioa; 

i ilenanes ; 

date; 

sound-speeda ; 
bathymetry ; 
wavefront ; 
paths ; 


all  scalar  parameters  in  initialization  file 
all  filenames  in  initialization  file 
date  and  time  of  ray  trace 
sound-speed  profiles 
bathymetry 

all  rays  at  receiver  range 
path  of  each  ray 


These  variables  and  their  contents  are  detailed  below  by  group. 


Initialization  Group 

The  initialization  group  contains  just  one  variable  named  inits.  It  has  a  size  of  I  by  25  and 
contains  the  version  number  and  all  of  the  scalar  parameters  (numerical  parameters  and  choice 
parameters)  contained  in  the  initialization  file.  Its  contents  are  detailed  below.  The  order  of  the 
contents  of  this  variable  always  match  the  order  in  the  prototype  initialization  file  created  when 
the  -Mksinit  command  line  flag  is  used. 


inits(l)  ss  ray  version  number 
inits(2)  s  source  range  (m) 
inits(3)  =  source  depth  (m) 
inits(4)  =  receiver  range  (m) 
inita(5)  =  receiver  depth  (m) 

inita(6)  —  first  angle  to  be  launched  at  source  (rad) 
inits(7)  =  last  angle  to  be  launched  at  source  (’ad) 
inits(8)  =  number  of  rays  to  be  launched  (integer) 
inits(9)  =  starting  range  for  saving  paths  (m) 
inits(lO)  =  ending  range  for  saving  paths  (m) 
inits(U)  =  fixed  step  size  along  a  path  (m) 
inits(12)  -  steps  per  path 


inits(13)  =  range  dependence  choice  (  0-3  ) 
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inits(14) 


inits(15) 


inits(16) 

inits(17) 

inits(18) 

inits(19) 

inits(20) 

inits(21) 

mlts(22) 

inits(23) 

inita(24) 

inits(25) 


0:  dc/da  range  dependence  (default) 

1 :  no  range  dependence  non* 

2:  dc/dz  range  dependence  grad.x 
3:  full  range  dependence  tall 
=  integration  routine  choice  (  0-3  ) 

0:  4th  order  runge  kutta  (default) 

1:  2nd  order  runge  kutta  rk_2 
2:  2nd-3rd  order  runge  kutta  rk.23 
3:  4th  order  runge  kutta  rk_4 
=  bathymetry  choice  (0-4) 

0:  no  bathymetry  (flat)  (default) 

1:  no  bathymetry  (flat)  none 
2:  soft  bottom  aott 
3:  absorbing  bottom  abaorbisg 
4:  reflecting  bottom  reflecting 
=  step  size  multiplier 
=  maximum  step  size  (m) 

—  minimum  step  size  (m) 

=  debias  factor  (  0  <  bias  factor  <  1) 
=  debias  iteration  (integer) 

=  bathymetry  smoothing  width  (m) 

=  bottom  depth  (m) 

=  radius  of  earth  (m) 

=  z  tolerance  (m) 

—  margins  (integer) 


Note:  Depth  scale  is  zero  at  the  surface  and  positive  downward  so 
that  all  depths  should  be  greater  than  or  equal  to  zero. 

Note:  Angles  are  defined  so  positive  angles  indicate  rays  moving  upward 
toward  the  surface. 

Note:  Some  of  these  parameters  may  change  before  Ray  runs  the  ray  trace, 
the  inits  variable  is  only  meant  to  reflect  what  the  user  has  input. 


Filenames  Group 

The  filenames  group  contains  the  names  of  all  of  the  filenames  specified  in  the  initialization 
file  plus  the  name  of  the  initialization  file.  Their  sizes  depend  on  the  number  of  characters  in  each 
filename.  They  are  all  MatLab  text  variables. 


variable  name  contents 


bathjfUe 

in  it -file 
mppJUe 
out -file 
prof-file 


bathymetry  filename 
initialization  filename 
MPP  input  filename 
output  filename 
profile  input  filename 
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Date  Group 


The  date  group  contains  two  variables.  One  is  a  text  variable  named  ron_date  which  contains 
the  date  and  time  that  the  ray  trace  started.  The  second  variable  is  named  runtime  and  contains 
the  number  of  seconds  that  elapsed  while  Ray  was  performing  the  ray  trace. 

Sound  Speed  Group 

The  sound  speed  group  contains  the  variables:  sspr,  idex,  ctabOOl  -  ctab###,  rec_c,  and  rtab. 
These  variables  contain  the  tables  that  Ray  uses  internally  for  modeling  the  speed  of  sound.  The 
sizes  of  many  of  these  variables  will  vary  depending  upon  the  data  in  the  input  hies.  Some  useful 
numbers  are: 

N„  The  number  of  sound-speed  profiles. 

Ni,at  The  number  of  bathymetry  sections  which  will  be  approximately 
twice  the  number  of  bathymetry  points  input). 

Nteg  The  number  of  segments  used  for  generating  sound  speed 
profiles  (usually  twice  to  three  times  the  maximum  number 
of  points  input  in  any  sound-speed  profile). 

JVj  the  number  of  meters  to  the  bottoa^lapth. 

The  variables  in  the  sound  speed  group  are  described  below. 
sspr(l  by  N„) 

The  range  (in  meters)  of  the  Nt,  sound-speed  profiles. 


idex  (Nt  by  1) 

A  table  of  indices  for  use  with  ctabOOl  -  etab###,  ree_c,  and  rtab.  For  every  meter 
of  depth  from  0  to  Nt  -  1,  idex(z  +  1)  tells  which  parameterization  of  the  ocean  to  use 
for  that  depth.  In  order  to  calculate  a  sound  speed  (or  a  step  size)  the  depth  in  meters 
is  used  as  an  index  into  this  array.  For  example,  to  find  a  sound  speed  at  a  depth  of  z 
meters  then  one  uses  ii  =  idex([z]  +1)  where  [z]  is  the  integer  part  of  z.  The  parameters  in 
etab###(ti  +•  1, :)  are  used  to  compute  the  sound  speed. 

ctabOOl  -  ctab###  (N,eg  by  4) 

Tables  used  to  construct  sound-speed  profiles.  The  numerical  suffix  will  vary  from  001  to 
N„.  They  are  ordered  by  increasing  range  from  the  source. 
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Column  1  is  the  depth  we  expand  about:  zq 
Column  2  is  the  sound  speed  at  zo:  eo 

Column  3  is  dc/dz  at  so:  ci 

Column  4  is  one  half  of  dPc/dz 2  at  Zo=  cj 
Sound  speeds  c(z)  are  calculated  as  follows: 

it  =  idex([r]  +  1)  +  1 
Zo  =  etab(ii,  1) 

Co  =  etab(tt,2) 
cj  =  ctab(u,  3) 
cj  =  etab(ii,4) 
dz  =  z  -  zo 

c(z)  =  c0  +  dzci  +  dz2  c  2 
tec.c  (N,eg  by  4) 

Table  used  to  calculate  the  sound  speed  at  the  receiver  range.  This  is  useful  since  the 
receiver  range  may  not  coincide  with  any  of  the  ranges  of  input  sound-speed  profiles.  The 
format  is  identical  to  that  of  ctab###  above. 

rtab  ( Nteg  by  5) 

Table  used  to  calculate  the  range  step  as  a  function  of  depth.  Its  use  is  similar  to  ctab 
described  above.  The  fifth  column  is  used  by  Ray  for  storing  auxiliary  width  information 
and  is  not  needed  for  computing  the  range  step  size. 


Bathymetry  Group 

The  bathymetry  group  contains  one  variable  btab  which  contains  the  table  of  values  Ray  uses 
for  computing  were  the  bottom  is. 

btab(JVfra(  by  5) 

This  table  is  used  to  determine  the  piecewise  parabolic  bottom.  The  ranges  in  column  1 
serve  a  dual  purpose.  They  are  used  as  the  range  about  which  the  parabolic  expansion  is 
made,  and  they  also  serve  to  mark  the  end  of  each  section  of  bathymetry. 


Column  1  is  the  last  range  of  each  section  r0 

Column  2  is  the  depth  at  ro  ho 

Column  3  is  db/dr  at  ro:  hi 

Column  4  is  one  half  of  d?b/dr 2  at  r0:  bj 


Column  5  is  the  minimum  depth  on  this  section  minz 
The  bottom  depth  h{r)  is  calculated  as  follows: 
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Find  the  entry  in  btab  such  that  btab(»  -  1, 1)  <  r  <  btab(t,  1) 
then  let: 
r0  =  btab(i,  1) 

6o  =  btab(i,2) 

6l=btab(i,3) 
b2  =  btab(»,4) 
dr  =  r  -  rO 

b(r)  ^  bo  +  dr  bi  +  dr 2  b2 

Note  that  the  number  of  rows  in  this  table  will  be  approximately  twice  the  number  of  input 
bathymetry  points.  If  there  is  no  bathymetry  specified  then  this  variable  will  not  exist. 

Wavefront  Group 

The  wavefront  group  contains  one  variable  called  wf  which  stands  for  “wavefront.”  Its  size  is 
Nr  by  20  where  Nr  is  the  number  of  rays  that  were  shot.  This  one  variable  contains  all  of  the 
output  from  the  ray  tracing.  Its  contents  are  described  below. 


wf(l)  =  angle  at  the  source  (rad) 

wf(2)  =  depth  at  the  source  (m) 

wf(3)  =  range  of  source  (m) 

wf(4)  =  angle  at  the  receiver  (rad) 

wf(5)  =  depth  at  the  receiver  (m) 

wf(6)  =  range  of  receiver  (m) 

wf(7)  =  travel  time  (s) 

wf(8)  =  number  of  top  turns 

wf(9)  =  number  of  bottom  turns 

wf(10)  =  number  of  top  bounces 

wf(U)  =  number  of  bottom  bounces 

wf](12)  =  total  number  of  turns 

wf(  13)  s  signature  flag** 

wfl(14)  =  upper  turning  depth  at  receiver  (m) 

wf(15)  =  lower  burning  depth  at  receiver  (m) 

wf(16)  =  eigen  miss  flag** 

wf(17)  s  arrival  type  flag** 

wf(18)  =  sound  speed  at  source  (m/s) 

wf(19)  =  closest  approach  to  bottom  (m) 

wf(20)  =  range  at  closest  approach  to  bottom  (m) 

Note:  **  indicates  variables  that  are  not  fully  implemented. 

The  upper  and  lower  turning  depths  are  depths  at  the  upper  and  lower  turning  points  that  would 
occur  if  the  ray  were  continued  past  the  receiver  range,  using  a  range  independent  sound-speed 
profile  identical  to  the  profile  at  the  receiver. 
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Paths  Group 

The  paths  group  contains  one  variable  named  paths  which  contains  all  the  path  information 
for  all  of  the  rays  that  were  traced.  The  number  of  columns  in  paths  depends  on  which  columns 
were  selected  in  the  initialization  hie.  If  all  of  the  columns  were  selected  then  paths  will  have  nine 
columns  and  the  will  be  in  the  following  order: 

paths(l)  =  range  (m) 

paths(2)  =  depth  (m) 

paths(3)  =  time  (s) 

paths(4)  =  angle  (rad) 

path8(5)  =  sound  speed  (m/s) 

paths(6)  =  dc/dz  (1/s) 

paths(7)  =  number  of  top  bounces 

paths(8)  =  number  of  bottom  bounces 

paths(9)  =  the  number  of  this  ray 

If  only  a  subset  of  the  possible  columns  are  selected  then  the  number  of  columns  in  paths  will 
change  but  the  order  will  be  as  shown  above.  The  number  of  each  ray  is  used  to  uniquely  identify 
the  rays.  The  first  ray  shot  is  numbered  one,  the  second  ray  is  numbered  two  and  so  on. 
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6.  Performance 


Range-Dependent  Environment 

Acoustic  data  from  a  3000  km  transmission  [Spieaberger  and  Metzger ,  1991]  and  the  Slice  89 
experiment  [ Duda  et  at.,  1992]  indicate  that  acoustic  wavefronts  are  smoother  than  the  predictions 
made  by  the  RDRYT  and  MPP  ray  tracing  codes.  The  smoothness  of  the  experimentally  measured 
wavefronts  is  taken  to  be  representative  of  ocean  acoustic  propagation  over  megameter  ranges  at 
frequencies  of  a  few  hundred  Hz. 

Figure  2  shows  the  arrival  time  predictions  for  runs  of  RDRYT  and  Ray  with  the  same  en¬ 
vironment.  Transmission  was  modeled  for  the  Pacific,  and  with  acoustic  source  at  1000m  depth, 
10.1°N  and  151°W.  The  receiver  was  2005.65  km  distant,  at  30°N  and  150°E,  at  a  depth  of  3000 
in.  Sound-speed  profiles  at  100-km  intervals  along  the  acoustic  path  calculated  from  the  Levitus 
temperature  and  salinity  database  for  summer  [ Levitus ,  1982].  This  database  has  profiles  at  one 
degree  intervals,  located  at  half-degree  positions  (e.g.  27.5°N,  68.5°W).  The  arrivals  predicted  by 
Ray  have  uniform  amplitudes  and  come  in  the  expected  groups  of  four.  The  RDRYT  arrivals  dis¬ 
play  the  groups  of  four  but  also  have  many  extra  arrivals.  The  bottom  half  of  Figure  2  shows  one 
of  the  groups  of  arrivals  with  an  expanded  time  scale.  It  shows  the  many  extra  arrivals  predicted 
by  RDRYT. 

The  nature  of  these  extra  arrivals  can  be  seen  more  clearly  in  Figure  3  which  shows  the  depth 
vs.  arrival-time  “timefront”  and  the  depth  at  receiver  range  vs.  launch  angle.  The  output  of  Ray  is 
smooth  and  evenly  spaced,  as  one  expects  for  a  sound  channel  which  is  slowly  varying  with  respect 
to  ray-loop  length  [Flatte,  1983].  The  output  from  RDRYT  has  many  zig-zags  and  an  uneven 
spacing  which  causes  the  extra  arrival  predictions  and  the  greater  variation  in  arrival  amplitude. 

The  depth  at  receiver  range  vs.  launch  angle  clearly  shows  the  contrast  between  the  smooth 
wavefront  generated  by  Ray  and  the  discontinuous  wavefront  created  by  RDRYT.  The  extrema  of 
these  curves  occur  at  caustics.  If  the  wavefront  is  smooth  in  the  neighborhoods  of  these  extrema 
then  it  should  be  possible  to  make  a  diffracted  extension  of  the  the  wavefront  in  the  shadow  zones 
of  these  caustics  ( Buchal  and  Keller,  I960].  It  seems  that  such  an  extension  would  be  meaningful 
with  the  results  from  Ray  but  not  with  the  results  from  RDRYT. 
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Figure  2.  Geometric  arrivals  for  a  2205.65  km  south  to  north  raytrace  in  the  tropical  Pacific  Ocean.  The 
source  depth  is  1000  m  tad  the  receiver  depth  is  3000  m.  The  lower  two  plots  are  expanded  versions  of  the 
second  complete  gtoup  of  four  arrivals. 


Flat  Layered  Environment 

In  a  flat,  range-independent,  layered  sound-speed  field  the  accuracy  of  numerical  ray-traces  can 
be  checked  against  analytic  solutions.  The  cosh(-r)  sound  channel  provides  an  interesting  analytical 
solution:  periodic  focusing  at  the  depth  of  minimum  sound-speed  (the  axis)  of  all  fully  refracted 
rays  emitted  from  the  axis  ( Tolstoy  and  Clay,  1966).  Calculations  were  made  using  the  model 
c(z)  =  a(cosh(6(z  -  *„)),  with  a  =  1480  m/s,  b  =  0.00006  nr1,  and  z0  -  2500  the  depth  of 
minimum  sound  speed.  The  foci  are  at  distances  nr/b  from  an  axial  source.  Figure  4  shows  a 
comparison  between  Ray  and  RDRYT  at  the  40th  focus,  which  is  at  a  range  of  approximately  2094 
km.  Values  of  c(z)  are  provided  to  Ray  at  100  m  intervals,  and  the  widths  u>,  we  all  50  n.  Note 
the  simultaneous  arrival  at  the  axis  (the  focal  points)  of  all  trajectories  traced  with  Ray,  agreeing 
with  the  analytic  result.  The  foci  were  missed  by  RDRYT,  seen  by  the  variable  depth  and  time  for 
rays  at  the  focal  distance. 


-0.10  -0.05  0.00  0.05 

launch  angle  (rod) 


0.10 


Figure  4.  Arrival  depths  and  times  for  rays  traced  in  a  cosh(z)  profile  for  exactly  40  focal  distances.  Each  ray 
traced  with  Ray  (solid  line)  passed  through  the  locus,  except  for  those  that  interacted  with  the  surface,  and 
they  all  arrived  at  the  same  time.  Trajectories  calculated  with  RDRYT  (dashed  line)  do  not  pass  through 
the  focus.  The  diamond  shows  the  analytic  arrival  time. 
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7.  Summary 

The  program  Ray  has  been  implemented  the  algorithms  of  Section  2  in  the  manner  described. 
The  Program  Input  and  Program  Output  sections  are  intended  to  provide  all  necessary  information 
to  use  the  program  properly. 

Ray  has  been  shown  to  achieve  its  goal,  the  calculation  of  smooth  acoustic  wavefront  propaga¬ 
tion  through  a  realistic  model  of  the  ocean. 
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APPENDIX  A.  Derivation  of  the  equations  of  motion 

The  equations  of  motion  for  a  ray  through  an  inhomogeneous  medium  can  be  derived,  by  a 
trivial  application  of  the  calculus  of  variations,  from  Fermat’s  principle  of  least  time  which  was 
expressed  by  Feynman  as  “  ...  a  ray  going  in  a  certain  particular  path  has  the  property  that  if  we 
make  a  small  change  ...  in  the  ray  in  any  manner  whatever,  ...  there  will  be  no  first  order  change  in 
the  time ."  [Feynman,  Leighton  and  Sands ,  1965]  Mathematically  this  idea  is  expressed  by  setting 
the  variation  of  the  time  to  zero, 


ST 


=  6  j  dt  = 


where  the  integral  on  the  right  hand  side  is  over  the  path  of  the  ray  with  fixed  limits  of  integration 
and  S  is  any  (differentiable)  variation  in  the  path  that  keeps  the  end  points  fixed.  In  order  to 
find  the  path  we  convert  the  integral  over  dt  to  an  integral  over  path  length  ds  using  the  index  of 
refraction  n  -  dtjds  -  1/c,  then 


ST 


=  6  J  nds. 


The  path  length  can  be  expressed  in  terms  of  cartesian  coordinates  x,  a s  ds  =  (dx.dxj)1/2  where  we 
have  used  the  Einstein  summation  convention  of  implicitly  summing  over  repeated  indices.  Take 
the  variation  of  this  equation  to  find 

Sds  =  dx,  Sdx ,  ds-1 
=  i,  Sdxi 
-  i;  Si,  ds 

where  a  dot  over  a  quantity  means  the  total  derivative  of  that  quantity  with  respect  to  s.  The 
variation  of  the  time  may  be  written  as 

ST  =  j  (S  n  +  niiSii)ds 

=  J  (6  Xidin  +  nx,6ii)ds. 

The  first  term  in  the  integrand  represents  the  change  in  time  due  to  the  change  in  the  index  of 
refraction  over  a  new  path.  The  second  term  represents  the  change  in  time  due  to  the  change  in 
the  length  of  the  path.  The  6  z,  dependence  of  the  integrand  is  eliminated  by  noting  that 

■^(niiSxi)  =  Sxi-^(nii)  +  nx,Sx,, 

and  integrating  the  second  term  in  the  integrand  by  parts.  The  result  is 

ST  =  J  Sx,  ^ di  n  -  d*- 

The  total  derivative  term  vanishes  due  to  the  assumption  that  the  variation  Sn  is  zero  at  the  limits 
of  integration.  The  equations  of  motion  can  be  obtained  by  demanding  that  this  expression  for  the 
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variation  of  the  time  be  zero  for  any  set  of  variations  6  This  is  possible  only  if  the  term  that  is 
multiplying  the  variation  vanishes  everywhere  along  the  path  of  integration.  Setting  this  term  to 
zero  gives  the  equations  of  motion  for  a  ray, 

This  may  be  rewritten  in  a  slightly  more  transparent  form, 


The  components  z,  form  a  unit  vector  pointing  along  the  direction  of  the  ray.  The  second  term  on 
the  right  hand  side  is  the  projection  of  Vn/n  in  the  x  direction.  Thus  x  is  equal  to  that  put  of 
Vn/n  which  is  perpendicular  to  the  path  of  the  ray. 

In  order  to  put  these  equation  in  a  form  suitable  for  (one-way)  numerical  calculations,  define 
a  horizontal  r  axis  and  a  vertical  x  axis  with  8  the  angle  of  the  ray  with  respect  to  the  horizontal. 
Then  x  =  f  cos  8  +  isintf,  and  the  equations  of  motion  become 

dd 

n(-f  sin  0  -(-  zcosB)  —  =  rdTn  +  zd,n  -  (f  c  os8  +  isin0)  (cos  9drn  +  sin  9dtn) 
ds 

The  r  component  of  this  equation  gives  an  expression  for  d9/ds , 

d9 

n-~  =  dr  n  sin  9  -  d.n  cos  9. 
ds 

Simple  geometry  gives  us  d9/ds  =  cos  8d9/dr  and  dn/n  =  -de/e.  Finally, 

-£l£. 

dr  c  e 

The  other  two  equations  of  motion  given  in  the  body  of  the  report  follow  from  the  definition  of  9 
and  simple  geometry. 


41 


APPENDIX  B.  Example  Initilisation  File 


/OMMOMMMMMMMO*  Ray  MMMMMMtMOMMOM/ 

/* 

•  Version  1.00  Friday,  loveaber  13,  1902,  3:33  pm 

•  Copyright  (C)  1992,  Woods  Hols  Oceanographic  Institution. 

•  Uss  -license  flag  for  ua«,  copying  and  distribution  conditions. 

*/ 


input  < 

/•  mpp.file  *  */ 

bath.file  *  "tsstl.bth"; 
prof.file  *  "tsstl.ssp"; 

>; 


output  { 

■at. file  •  "test  1  .Mt" ; 
initialization; 

filenames; 

date; 

sound.speeds ; 
bathymetry; 
paths ; 
wavefront ; 

/*  everything;  •/ 

/•  environment  .only ;  ♦/ 


source  { 

range 

depth 

>; 

receiver  { 
range 
depth 

>; 


0  (km); 
1000  (m); 


2000  (km); 
0  (m); 


angles  { 

first 

last 

number 
/•  specific 

>; 


>  IS  (degrees) ; 

■  -IS  (degrees); 

■  800; 

■  {}  (degrees);  •/ 


paths  { 

■in. range 
nax.range 
fixed. dr 
/*  steps. per 
columns  { 

range; 


0  (km); 
2000  (km); 
600  (■); 
0;  */ 
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depth; 

time; 

•ngls; 

■paad; 

grad; 

top. bounces; 
bo t .bounces; 
ray .nunber ; 
/•  everything; 

>; 

include  bounces ; 

}; 


nodel  { 

range. depend 

Integration 

bathynetry 


grad.z;  /•  none  grad_z  full  •/ 
rk.4;  /•  rk_2  rk_23  rk_4  •/ 

nose;  /•  none  soft  absorbing  reflecting  */ 


range. step  { 

Multiplier  • 

aaz  ■ 

Min  ■ 

>; 


0.5; 

200  (■); 
5  (■); 


>; 


dsbias  { 

factor  ■ 
iteration  ■ 

>; 

/•  prof .smoothing 
bath.snoothing 
bottoa.depth  ■ 
earth.radius  « 
z.toltrance  ■ 
Margins  ■ 


i; 

10; 

<>  (■);  •/ 

10  (kM); 
6000  (n); 
6378.137  (kM); 
le-06  <m); 
100; 
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APPENDIX  C.  Input  File  Formats 


Profile  File  Format 

A  profile  file  contains  the  information  necessary  for  Ray  to  compute  the  sound  speed  at  all 
ranges  between  the  source  and  receiver.  It  is  organized  as  a  sequence  of  tables  of  sound  speed 
(m/s)  as  a  function  of  depth  (m)  at  increasing  range  (km).  For  each  sound-speed  profile  there  is  a 
header  line  containing  1)  the  range  of  the  profile  in  km,  2)  the  number  of  points  in  the  profile,  3) 
the  number  “0”  Directly  after  the  header  is  a  list  of  the  depths  (m)  and  sound  speeds  (m/s),  with 
one  pair  of  numbers  on  each  line.  After  the  listing  of  the  profile  points  there  is  a  single  line  with 
the  a  “0”  on  it.  This  pattern  of  header  line  followed  by  profile  points  followed  by  a  “0,”  is  repeated 
until  the  end  of  the  file  which  is  indicated  by  the  string  “END”  following  the  single  “0.”  Here  is  a 


short  example. 

0.0 

33 

0.0 

1517.4268 

10.0 

1617.6736 

20.0 

1617.8807 

30.0 

1618.0923 

so.o 

1618.3930 

75.0 

1616.7801 

100.0 

1611.6797 

126.0 

1608.6289 

160.0 

1607.2884 

200.0 

1604.7734 

260.0 

1602.6900 

300.0 

1600.0790 

400.0 

1496.4920 

500.0 

1489.6708 

600.0 

1404.1949 

700.0 

1480.7405 

900.0 

1479.4832 

900.0 

1479.4336 

1000.0 

1479.7936 

1100.0 

1480.6996 

1200.0 

1481.5168 

1300.0 

1482.6431 

1400.0 

1483.5417 

1600.0 

1484.6309 

1760.0 

1487.3678 

2000.0 

1490.5309 
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2600.0 

3000.0 

3600.0 

4000.0 

4600.0 

6000.0 

6600.0 

0 

220.0 

0.0 

10.0 

20.0 

30.0 

60.0 

76.0 

100.0 

126.0 

160.0 

200.0 

260.0 

300.0 

400.0 

600.0 

600.0 

700.0 

800.0 

900.0 

1000.0 

1100.0 

1200.0 

1300.0 

1400.0 

1600.0 

1760.0 

2000.0 

2500.0 

3000.0 

3600.0 

4000.0 

4500.0 


1497.8366 
1605.9020 
1614.4027 
1523.1536 
1532.1662 
1641.2716 
1660.5268 

32  0 

1516.7812 
1616.0229 
1616.1872 
1616.3413 
1616.6648 
1613.7440 
1608.9766 
1606.3369 
1601.5737 
1497.8691 

1496.8367 
1493.9893 
1488.9931 
1483.0016 
1479.2471 
1478.6666 
1478.8672 
1479.3947 
1480.1660 
1480.9617 
1481.9212 
1482.9679 
1483.9375 
1484.9661 
1487.4470 
1490.6247 
1497.8800 
1606.8644 
1614.3496 

1623.1794 

1632 . 1795 
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$000.0 

OEID 


1641.3340 


Bathymetry  File  Format 

The  bathymetry  file  specifies  the  depth  of  the  bottom  as  a  function  of  range.  The  format 
consists  of  two  columns  of  numbers  specifying  the  range  (km)  and  the  depth  (m).  The  number  of 
bathymetry  points  is  equal  to  the  number  of  lines  in  the  file.  Here  is  a  short  example. 


0.00 

1000.00 

10.00 

1200.00 

20.00 

1400.00 

30.00 

1600.00 

40.00 

1800.00 

60.00 

2000.00 

00.00 

2200.00 

70.00 

2400.00 

80.00 

2600.00 

00.00 

2800.00 

100.00 

3000.00 

110.00 

3200.00 

120.00 

3400.00 

130.00 

3000.00 

140.00 

3800.00 

160.00 

4000.00 

180.00 

4200.00 

170.00 

4400.00 

180.00 

4600.00 

190.00 

4800.00 

200.00 

6000.00 

210.00 

6200.00 

220.00 

6400.00 
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